VIB/IR analysis workflow using MatSciToolkit#
This tutorial will guide you through the use of the MatSciToolkit package to calculate the vibrational frequencies of a molecule or surface. The MatSciToolkit package is a python package the provides a set of tools to automate the workflow for the vibrational and IR analysis. The package is built on top of the Atomic Simulation Environment (ASE) and simplifies the workflow while using the Quantum Espresso code as the calculator.
A quick example of the data we get is the following
Show code cell content
from ase.build import molecule
from ase.calculators.emt import EMT
from ase.optimize import BFGS
from ase.vibrations import Vibrations
from ase.visualize.plot import plot_atoms
import matplotlib.pyplot as plt
import matplotlib.animation as animation
h2o = molecule('H2O', cell=(10,10,10))
h2o.center()
h2o.calc = EMT()
BFGS(h2o).run(fmax=0.01)
Show code cell output
Step Time Energy fmax
BFGS: 0 16:24:28 2.619811 7.7384
BFGS: 1 16:24:28 1.912906 1.3454
BFGS: 2 16:24:28 1.882033 0.4035
BFGS: 3 16:24:28 1.879275 0.0317
BFGS: 4 16:24:28 1.879253 0.0096
True
Show code cell source
import plotly.graph_objects as go
from ase.data.colors import jmol_colors
from ase.data import atomic_numbers, covalent_radii
import numpy as np
def ms(x, y, z, radius, resolution=20):
"""Return the coordinates for plotting a sphere centered at (x,y,z)"""
u, v = np.mgrid[0:2*np.pi:resolution*2j, 0:np.pi:resolution*1j]
X = radius * np.cos(u)*np.sin(v) + x
Y = radius * np.sin(u)*np.sin(v) + y
Z = radius * np.cos(v) + z
return (X, Y, Z)
pos = h2o.get_positions()
numbers = h2o.get_atomic_numbers()
colors = [jmol_colors[number] for number in numbers]
cs = np.array(colors) * 255
cs = [f"rgb({int(color[0])}, {int(color[1])}, {int(color[2])})" for color in cs]
scale = 1
radii = covalent_radii[numbers]
radii = [(scale * r)**1 for r in radii]
data = []
for i in range(len(pos)):
(x_pns_surface, y_pns_surface, z_pns_surface) = ms(pos[i,0], pos[i,1], pos[i,2], radii[i])
data.append(go.Surface(x=x_pns_surface, y=y_pns_surface, z=z_pns_surface, opacity=1, showscale=False, colorscale=[cs[i], cs[i]]))
fig = go.Figure(data=data)
fig.update_layout(
title='Water molecule structure',
scene=dict(
xaxis=dict(title='X', range=[0, 10]),
yaxis=dict(title='Y', range=[0, 10]),
zaxis=dict(title='Z', range=[0, 10]),
aspectmode='cube'
),
width=1000, # Set the width of the figure
height=500, # Set the height of the figure
margin=dict(l=0, r=0, t=30, b=0), # Set margin to 0 for a tight layout
template='plotly_dark',
)
fig.show()
We’ll get the vibrations or infrared data#
vib = Vibrations(h2o)
vib.run()
vib.summary()
---------------------
# meV cm^-1
---------------------
0 6.3i 51.0i
1 5.9i 47.6i
2 0.0i 0.3i
3 0.0i 0.1i
4 0.1 0.4
5 5.4 43.8
6 32.1 258.9
7 296.7 2392.7
8 387.4 3124.9
---------------------
Zero-point energy: 0.361 eV
Show code cell source
vibdata = vib.get_vibrations()
To visualize, the mode 6, 7, and 8 is:#
Show code cell source
%matplotlib notebook
from IPython.display import HTML
plt.ioff()
fig, (ax1, ax2, ax3) = plt.subplots(3, 2, figsize=(8, 8), dpi=150)
axs = [ax1, ax2, ax3]
itervib6 = list(vibdata.iter_animated_mode(6))
itervib7 = list(vibdata.iter_animated_mode(7))
itervib8 = list(vibdata.iter_animated_mode(8))
itervibs = [itervib6, itervib7, itervib8]
fig.tight_layout()
def animate(i):
for ax in axs:
ax[0].cla()
ax[1].cla()
for k, (ax, itervib) in enumerate(zip(axs, itervibs)):
plot_atoms(itervib[i], ax=ax[0], radii=0.9, rotation=(''))
plot_atoms(itervib[i], ax=ax[1], radii=0.9, rotation=('90y'))
mode = {0: "Mode 6", 1: "Mode 7", 2: "Mode 8"}
ax[0].set_title(mode[k])
ax[1].set_title(mode[k])
# if k == 2:
ax[0].set_xlabel("X-axis")
ax[1].set_xlabel("Z-axis")
ax[0].set_ylabel("Y-axis")
ax[1].set_ylabel("Y-axis")
if i == 0:
fig.tight_layout()
ani = animation.FuncAnimation(fig, animate, frames=len(itervib6), interval=25)
# from IPython.display import HTML, Javascript
HTML(ani.to_jshtml())
# HTML(ani.to_html5_video())